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1. Introduction 


After many years of development of the hardware, the airborne gravity gradiometry has 
reached the operational stage. The test flight was taken in the Texas-Oklahomu area and the test 
results were published (Brzezowski, et al., 1988). The test area is very smooth, so the 
topographic effect was neglected. In the future, the airborne gravity gradiometry will be used for 
the rough mountain area. The effect of the topography has to be taken into account in more rugged 
topographic areas. For many years this problem has been studied by different authors (Chinnery, 
1961; Dorman et a!., 1974; Hammer, 1976; Tziavos, et al., 1988). All studies had a basic idea - 
they intended to eliminate the effect of the topography by removing the mass above the geoid. The 
gradients of the attractions of the mass above the geoid were subtracted from the aerial gradient 
data, and the gravity disturbances could be determined on the geoid by processing the reduced 
aerial gradient data. 

If the gravity disturbance is determined on the earth's surface, other methods can be used. 
One of the methods was suggested by iekcli (1987). He used a surface integral to determine the 
disturbing potential on the earth's surface and avoided using the topographic reduction. 
Theoretically, this method is perfect but it is difficult to realize in practice, because the inclination 
of the topography, which is ill-defined, is needed at every computation point. 

An alternative solution (ibid., p.239) which is difficult but simply defined is the use of the 
analytical continuation method. Assume that the derivatives of the disturbing potential T, such as 
T/. T//, T 777 , ... can be well determined at a mean plane through the topographic surface. By 
analytical downward continuation, the gravity disturbance can be determined on the topographic 
surface by using Taylor's series. 

It was shown (Schwarz, 1979; Rummei, et al., 1979; Neyman, 1985; Ilk, 1988) that 
downward continuation is an improperly posed problem. An improperly posed problem may have 
a solution but it does not depend on the data continuously. A small error in data, e.g. a random 
measurement error, can cause a significant deviation in the solution. It is expected that the second 

derivatives of the disturbing potential T, such as T 7-Z , T /x , T xy .are rough at the Bight altitude. 

The problem is, how can we downward continue these rough functions to a mean level? 
Furthermore, how can we absorb the useful information from such data to determine even higher 
derivatives of the disturbing potential on the mean level? Sometimes it looks like it is impossible, 
but if the gradient data is accurate and in good distribution, the reasonable results can always be 
expected. In order to avoid the instability of the computations and get a reasonable smooth 
solution, there are different methods that can be used, e g., least squares collocation, 
regularization, or smoothing (filtering). These methods have the same property: they filter out the 
high frequency ot the data and make the results stable and smooth. We will show that the three 
methods arc identical under some conditions . 

The goal of this study is a find methods for the determination of the gravity disturbance on 
the topographic surface by processing the aenal gradient data. The numerical computation will be 
carried out to gam an idea about the use of the methods. 

Because the gradient data can be obtained at regular gnd points, the very efficient numerical 
computation method Fast Fourier transformation (FFT), is used. We will study the problem in 
the spectral domain and use FFT in the nurnenca! computations. 


7 solution of Analytical 'Jownward Continuation for the Airborne Gravity Gradiometry 

This section presents the formulas of the analytical downward continuation for the airborne 

eras its gradiometry 
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Because the airbonre gravity-gradiometry is taken in a local area, the flat-earth approximation 
is suitable for the processing the aerial gravity gradient data (Jekeli, 1985). 

At first we consider the analytical downward continuation of the aerial gravity-gradient data to 
the mean elevation level. The geometry of the airborne gravity-gradiometry is drawn in Figure 1. 



Figure 1. Geometry of the Airborne Gradiomecry 


Assume that the Runge's theorem (Moritz, 1980, p. 67) is valid also for the plane approximation, 
one can say that there is a function which is harmonic on and outside the mean elevation level and 
this function approximates the disturbing potential on and outside the earth's surface as good as we 
wish. We assume this function can be approximated by the analytically downward continuing the 
disturbing potential of the earth from outside the earth's surface to the mean elevation level. 

Therefore we assume that the disturbing potential T and its derivatives, such as Tj, Tjj, i, j = 
1. 2. 3 corresponding to the subscripts x, y, z respectively, are analytically downward continued 
inside the earth and are harmonic above the mean elevation level. 

The Poisson's integral gives the relationship between a solid harmonic function and its values 
at the mean elevation level (cf. Heiskannen and Moritz, 1967, p. 239): 


Tf j Z H f f T\x, y) 

T (Vy P ’G=—JJ -^dxdy 

2n J J i r 


(1) 


vshere 1 = [fx-x p ) 2 + (y-yp) 2 + zfi| 1/2 , is the distance between the current point on the mean 
elevation level t and the point P at the flight level; z\\ is the height of the flight level above the mean 
elevation level. Eq. (1) is valid for any harmonic function. For the derivatives of the disturbing 
potential T,, T,j we have: 
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The second derivatives of the disturbing potential Tu are given on the flight level, the 
derivatives of the disturbing potential T, such as Tj, T 1Z , Tj TL and even higher terms will be 
determined on the mean elevation level. We consider this process in two steps: first the 
components of the gravity disturbance Tj are computed at the flight level by processing the second 
derivatives of the disturbing potential T,j The formulas can be found in (Jekeli, 1985) 
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( 6 ) 


(7) 


( 8 ) 


where If, = (x-x'j + (y-y')2 and u is the azimuth of the point (x', y') with respect to the point 

I x, y ). 

The second step is to downward continue the derivatives of the disturbing potential Tj, T 1Z , 
T t// to the mean level by using the Poisson s integral. If Tj, T 1Z , T VJJ are determined on the mean 
level, then we can get the gravity disturbance on the topographic surface by using Taylor's series: 


T, |Q) = T,(p'! + Ah r ' 


r/r. 


f)h 


h = hm 1 


-(Ah)' 
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a t. 


dh 
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= T, (p'l v Ah T 1/ (P'|+ - (Ah)' T l7/ (p'j 
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i = 1.2,3 


(9) 
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w here Ah = QP' = hg - h m , the height of the topography referenced to the mean elevation level. 


The role of the mean elevation level is like the point level in the analytical downward 
continuation solution of the Molodensky's problem (Montz, 1980). Instead of the point level in 
the solution of the Molodensky's problem the mean elevation level is being used, so that Tj, Ah 
T|/, 1/2 (Ah)- T 17Z should have a similar magnitude and property as Ag, gi and g2 of the solutio.. 
of the Molodensky's problem. For more details about the numerical properties of gi, g 2 and 
higher terms see (/Wang, 1987). 

In the following discussion we consider how the gravity disturbance and its derivatives can 
best be determined on the mean elevation level. First we consider the use of least squares 
collocation to process the aerial gravity-gradient data. 


.V Processing the Aerial Gravity-Gradient Data by Using Least Squares Collocation in a 
Continuous Case 

If the data are dense and regularly distributed as in the case of airborne gradiometry, least 
squares collocation can be considered in a continuous case. The advantages are that the problem 
c m be solved in the frequency domain easily and the efficient numerical computation method - Fast 
Fourier transform can be used. 

Generally, see consider the operator equation 

g' = Af g° e G, f e F (10) 

a here A: F —> G, a linear operator which maps the normed space F into the normed space G. In 
airborne gradiometry A can be any integral or differential operator. A specific example is: A is the 
upward continue operator defined by the Poisson's integral (formula (3)); f is the second derivative 
of the disturbing potential T,j at the mean level T and g° is T,j at the flight level. Now we have the 
second derivatives of the disturbing potential T,j at the flight level. We want to determine T,j on 
the mean level I. This inverse problem may be properly posed if T,j at the flight level is smooth 
enough and errorless. 


In practice such an inverse problem is improperly defined, because we always have the errors 
m the data. It means, instead of the original function g°, we have in practice 

(II) 

•where r is the measurement error 

IN en though the inverse of the operator A ' exists, the solution 

f* - A' 1 g (12) 

vim he qu.ie different from t which we are trying to determine. Now we want to find the methods 
to overcome (has difficulty 

It we have (he previous information about the statistics of the data and the measurement error, 
the method of least squares collocation can he used to obtain a solution. This technique has 
become a standard computational method for the inverse problem in physical gecxlesy. 








In eq. (11 i vu: assume that the function g° is centered: 
M |g 0 } =0 


The operator N1 is defined as 


M (g c / = lim 

T—» 


i r T r T 

j l g° (x, vj dx dy 
4T“ J T J i 


( 13 ) 


The function g° is deterministic and the error £ is considered as a stochastic process, so the 
measurement g is a stochastic process. 

The variance and the covanance funcnon are defined as 

Cff|P, Q) = M {f (P) flQ)} 

(14) 

C rt ,(P,Q) = M{f(P|g(Q} 

P. Q are the points on the reference plane. 

We assume the function g° and the error £ are independent: 

M {eIP) g°!Qj = M (gIP)eIQ)/ =0 

M (e!P) gjQ)) - M (e |P) £ |Q)j --- C nn (P,Q| (l5) 

where C nn is the covariance of the error £. 

We consider a process for the best estimation of the function f: 

= (16) 

where f is the best estimate of the function f, H is the estimation operator which makes the mean 
square estimation error e the smallest: 


M fe 2 ) - M {(f t) 


min. 


(17) 


r.q ; 17; is equivalent to 

(',c (r - 0, H i - min. (18) 

where C cc is the covanance of the estimation error, and it is a function of the estimator H; r is the 
!; ounce between the points P and Q. 
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Eq. (18) can be viewed as an extreme value problem: To find an estimator H which makes 
C\v (r = 0, H) the smallest. 

Now we consider how to solve this minimum problem in die frequency domain. 

The two dimensional Fourier transformation and its inverse are defined by 


FlfU.yD^JJ'fS.ylO^'-'dxdy 

F ‘ Iglu.v)} = j j g ixi .X |e "' 1 ' i,u ' >vl du dv 

w here j =■/-!; F and F 1 denote die Fourier transformation and its inverse respectively. 
We denote the Fourier transformation of the function f by oif 

to- = F {ft v i 1 
and assume dial 

F | Af} = 0,\ cof 


(19) 


( 20 ) 


( 21 ) 


( 22 ) 


F lilt! ■- Cat (Of 

w here 0,.\. Oh are railed the spectra of the operators A and H respectively. 

For the Fourier transformation of the covariances we have (ct. Schwarz, et al., 1989): 


R „! U.s j 


= lm, —r(cu fM ;) 

T 4T 

(23) 

k , j u.s 

, = F { C .8 (X - V| ] 

[ = ! ,m ») 

T * » 4T 



hm —(a 
" •- 4 1" 


(24) 

R U. V 

1 : !C „ -x. yi, 

1 - o..\ Rif 



K*, 


(25) 

Rjj 

' 0,\ 0a kll 


(26) 
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where the symbol " * " denotes the conjugate of a complex function and Rff, Rfg, R g f, Rg C are the 
Fourier transformation of the covariances. Sometimes they are called the power spectral density 
function. 

If we denote the power spectral density of the estimation error by R^ (u, v), then eq. (18) can 
be written in the form. 


C,.Jr = 0, HI 


■//; 


Rec(u. V. 


(jjjJ du dv = min 


(27) 


then the extreme value problem becomes: 

To find the spectrum of the operator H which makes the estimation error the smallest. 

If there is a procedure which minimizes the power spectral density function of the estimation 
error everywhere in the frequency domain: 


R _(u. v, OpJ = min.. 


l!,V€ 


- oo foo | 


(28) 


the word "everywhere' means the minimum • aiues of Rc c for every frequency u, v, then the 
extreme vu'ue problem ( 18j can be replaced by (28) in frequency domain. The minimum condition 
Iv was used by Rendat (1980) for minimizing the power spectral density function of the 
estimation error. 


Note that the R cc is non-negative and compare eq. (28) with eq. (27), one can find that eq. 
I'i can be obtained by using eq. (28). If the power spectral density function R cc is minimized 
everywhere, its integration C cc (r - 0 H) is also smallest. 

The estimation error covariance is given by 

f:.= M 'e IP) e |Qf) = M {[f |P| - ?(P)] [f|Q - f !Ql]} 

= C,f - C, r C ff * Cn ^ 29 ) 

Using eqs. (23/ - (26) and /29) we get the power spectral density function of the estimation 
err< >r e: 


. I 


Rff Rff R ff + l^A ^A R ff + ^nn. 


(30) 


1 nder (he minimum condition (28), the spectra Oh. Oh have to satisfy the following conditions: 


') 

— |R ,J (| 
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c)0 H 


From 1 30") and (31) we get: 


R ff 

0 H - - 

0 A 0 A Rff + R nn 


( 32 ) 


0 H = 


0 A R ( 


o A o A Rff + R nn 

lie best estimate of the function f is given by 


* = ' : '!«„%) 


* ^ ! 
1 0 \ Q \ R ff * R i;n I 


chere to*. is the Founer transformation of the measurement (data) g. 
The estimation error is given by (of. eqs. (27) and (30)): 

M |e-) = C oc dll 


(33) 


(34) 


If 


1 0 A 0 H j Rff + R, 


dudv 


(35) 


If the pre information of the statistics of the function f and the error e, or the power spectral 
tensity function Rff and R nn . are known, then we can use formulas (32) and (34) to get the best 
estimation <>t the (unction f from the data g. The estimation error can he computed by using eq. 

i ' 5 . 


Now we consider a more general case. Wc have a heterogeneous data set related to the 

tun*.is>n f 
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gl = Ai f + E| 
g2 = A 2 f + t'2 


gn = A n f + E n . (36) 

A solution which is the best approximation of the function f is to be determined by processing 
this heterogeneous data. An example for eq. (36) is the processing of aerial gravity-gradient data, 
gi, i = 1, 2. ... 6 is the second derivative of the disturbing potential, and f can be the disturbing 
potential or any of us derivatives. Of course, if we have any other data related to the function f, we 
can always put it in the form as eq. (36). 

We assume the best estimate of the function f is given by: 


: = H : g, + H : g, + ... + H n g n = £H,g, 


i=i 


The estimation error is 


it 

e = f - f = f - H, 


(37) 


(38) 


The estimation error covariance has the same form as eq. (29): 

C ee =<VCf f .C f r + Cn ^ 29 ') 


Note that 



= X R II 

i=l 


R,f= F |C 
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n n * / * 

= Z Z^H^H l^i R ff +R i 

«=1 J=> 

ss here 

R _/<V° h 

1 ° '* J (40) 

Here we have assumed the measure errors e,, i = 1,2, ... n are independent of each other, Rjj is 
the power spectral density funedon of the measurement error £,; 0, is the spectrum of the operator 
•V The spectral density function of the estimation error is: 

n n • « 

R C e= R fr Z^tA R rr Z$h$j R tf + 

.= t ' j=t 


+ Z I VMM, R ff R ij 

1=1 J= l 


(41) 


The power spectral density function has the extreme value when: 

dR ,. 

—4i = 0 


dO, 


3R_e. 

dO, 


= 0 


Using eqs. (41) and (42) we get 


; ^ it • | 

Zo, k it + ZZ°h o .°, R ff +R u / =0 


'=!J=» 


n n * 


- Z a R ff + ZZ aiIaa R rf+ R 

>=ij=i 


i=i 


J =0 


Kq (43) can be rewritten as 
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0, R rr + Z^hK^j R ff + R 


i=i 


= 0 


(42) 


(43) 


(43*) 
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i = 1 


o, Rfr+ZonK^j R if + R ij, 

J=I 


= 0 


Eq. (43*) implies that 

* n I * > 

~ <+>_, Rff + X R ff + R J = 0 


i=l 


R ff + X OhKOj R ff + R J- o 

J=1 


(44) 


(45) 


Equations (44) and (45) are conjugate to each other. They are indeterminate linear equations. 
There are infinite solutions for the "variables" 0 h, and 0Hj- One of the solutions of eq. (44) can be 
written as: 


0, Rff 

0 H =- when R,, = R 22 = ... R nn . 

n • 

X Oj0, Rff + R,, 

, = 1 ' (46) 


This is a special case in which all measurement errors have the same magnitude and property 
(same power spectral density function). One can expect that it can happen sometimes, e.g., in the 
airborne gravity-gradiometrv, the components of the gradient of the gravity are measured with 
different accuracy. 

In the last case the solution for the indeterminate equation (44) can be 


°H = 


0! R ff R ll 


°, R if R jj +1 


(47) 


If the power spectral density function of the measurement error is "white" noise, we have: 

R,i = constant (48) 

We then can define a spectral weight oby 
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(49) 


W p 



so that equation (47) becomes: 


0 


H, - 


<\ Rffttp 


n . j 

0 J Rff«p+ 1 
j=i 


(50) 


Obviously, eq. (50) can be considered as a weighted least squares collocation solution with weight 
If a data set is measured with low accuracy, based on eqs. (49) and (50), the spectral weight 

becomes smaller and so the 4 >h,- This data set is weighted and has less contribution to the 
results. 


The definition of the spectral weight o)p can also be expanded to a more general case in which 
the R cc is not restricted to be a constant. Assume that the power spectral density function of the 
measurement errors is not only the "white" noise and let eq. (49) still be valid, eq. (47) can be 
considered as a weighted least squares collocation solution with weight (i)p, too. The only 
difference from the "white" noise case is that the spectral weight C 0 p has different values to different 
frequencies. 


We now consider only the case in which all data are measured with the same accuracy. The 
best estimation of the function f is given by 


f = F 



= F 




X +R n J 


(51) 


■a here ojg i is the Fourier transform of the data g,. 

The estimation error variance is (cf. eqs. (27) and (41)): 
M (e~) =Qc (0) 
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Here we have used: 


n n * 

IXW‘1 

i=lj=! 



( 53 ) 


In an ideal case the data are errorless, the spectrum of the best estimation operator H, has the 

form: 


n * 

I 0,°, 

1 = 1 


(54) 


then the function f can be exactly recovered without using any pre-information of the statistics of 
the function f. 

Because eqs. (44), and (45) are indeterminate linear equations and pose infinite solutions, we 
can find another solutions of the operator Hj which satisfies eqs. (44) and (45). An alternative 
method to solve this problem is discussed in (Bendat, et al., 1980, Chapter 10). 

The above formulas can be used for any observed quantities which are related to the disturbing 
potential, e.g., the gravity anomalies; deflections of the vertical; geoid undulation etc. If we have 
such heterogeneous data, the above formulas can be used to determine the needed quantities. 

The weakness of this method is that the data have to be regularly distributed and should not be 
so sparse that the interested information is lost. 


4. Regularization 

In Section 3 we have considered the solution of the improperly posed problem by using least 
squares collocation. Now we study the problem from another starting point. We will consider the 
solution to be stable and smooth. 

In least squares collocation the statistics of the determined function and the measurement error 
have to be known. In practice, they are never exactly known and are always assumed. If the 
solutions are sensitive to the statistical model of the measurement error or of the function being 
estimated, or the statistical model is not properly assumed, the results may not be good or not 
stable. In this case one can consider the use of a regularization method. 

The regularization method has been used in many technical and scientific areas (Nashed, 
1974) The use of the regularization in physical geodesy can be found in (Schwarz, 1979; 
Neyman, 1985; Ilk, 1988). This method is flexible in obtaining a stable and smooth solution 
because we have the chance to choose the regularization parameter and the regularization function 
arbitrarily. 


- 13 - 








There are different regularization methods and different ways to regularize an improperly 
posed problem; for more detailed see Nashed (1974). Here we are interested in this problem: Find 
a solution where its n 1 * 1 derivatives are smooth and it is the best approximation of the original 
solution. 

For the improperly posed operator equation 


g = Af 


( 10 ’) 


f e F, g e G 


Let F, G and Z be Hilbert spaces and A: F-»G be an operator mapping the Hilbert space F into G; 
L m : F—>Z be an operator mapping the Hilbert space F into Z. We consider the minimization 
problem: find a function f a to minimize the functional 



w here the norm of a function f is defined 

'irMf.o (56) 


Eq. (56) is from the definition of the Hilbert space (Bachman, et al., 1966, p. 141). The inner 
product (f, 0 can be defined for our purpose as: 


(f.f) =lim — ~ f f ff’dxdy 

T A T * J T J T 


T— 4T -T -T 


(57) 


where f* is the conjugate of the function f. 


Notice that 


Af - g = e 


(58) 


where e is the measurement error, therefore eq. (55) is equal to 
j(f, g.a, L m ) = ||e|| c +a j|L m f|j ? 


(59) 


Let a = 0, then the minimization problem (55) becomes: find a function f a to minimize the 
functional 


( 60 , 

Eq. (60) is a classic least-squares minimization problem. The physical meaning of the last term in 
(59; is clear. If L m is chosen as a differential operator up to m order, then the minimization 
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problem (55) means: to find a function f a to minimize not only the measurement error, but also the 
functional I tL m f I!. The last term I |L m f 11 makes the function f a smooth and stable and it makes 
the difference between the classic least-squares minimization problem and the regularization 
problem. 

The solution of (55) has been given by (Nashed, 1974, chapter 4): 

(a*A + cx L^L m )f a =A*g 

where A* is the adjoint operator and is defined by (Bachman, et al., 1966, p. 16): 

|Ax, y) = (x. A* y), x.yeF 
If the inverse of the operator 

(a* A ^ a ClJ 
exists, then we have the solution 

f rl — (a A + ct L m L m ] A g 

Based on the Lemma 2.2 (ibid. p. 25) we get the spectrum of the adjoint operator A* 


(62) 


(63) 


°a'-°a (64) 

That means that the spectrum of the adjoint operator A* is equal to the conjugate of the spectrum of 
the operator A. 

Applying the Fourier transformation to eq. (61), we get the spectrum of the regularization 
solution f a : 


F ( f„! = 0) f = 


i ,l l 


$.01 
A g 


0 A 0 A + ct 4>[^ 0 L 


(65) 


We denote (61) by an operator equation 




( 66 ) 


with 
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-1 


Hn = A A + a L m L 


ml 


Hq was viewed as a regularization operator by Nashed (1974). 

In the spectral domain eq. (66) can be written in the form (cf. eq. (65)): 




(67) 


( 68 ) 


with 


• - ; 0 2 m 

p A 0 A +a o L- o Ll , 

The regularization error is defined as 
e r = f - fa ■ 

and the error covariance is 


(69) 


(70) 


C[. c = M (e r {P)e r {Q)J 

Cff C fn (' 6- ffa + ^-fafa 

The power spectral density function of the regularization error is given by: 




-> 



where R nn is the power spectral density function of the measurement error. 


(71) 


(72) 


From eqs. (61) and (69) we can see that the solution f a is dependent on the regularization 
parameter a and on the choices of the operator L m . No statistical model of the function f and the 
measurement error are needed. But if the regularization error is to be determined, the statistics of 
the measurement error and of the function being estimated, are always needed (cf. eq. (72)). 

The mean square of the regularization error is given by: 




R r |u, vjdudv 
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( 73 ) 


-If 


1 l 2 

1 l 2 

l 1 ' °H 0 < > a) 

Rff+JOnJ Rnn 


dudv 


Q c (0, a, L n1 ) is a function of the regularization parameters cc and the operator Lm- 

The regularization of the improperly posed problem can be extended to the heterogeneous 
data. Denote the equations (36) with the vectorial form: 


g = A f 


(74) 


where the underbar denotes the vector, g is a vector function with the components g;, i = 1,2, 
...n, and A is a vector operator with components A,, i = 1, 2, ... n. The minimization problem is: 
Find a function f (1 , to minimize the functional 


J(f, g, a, L m ) = 11 A f-gil G + a 2 ||L m fllz 

The norm of a vector in the Hilbert space G is defined by the inner product: 

|!g|r = (g.g) = (Af.AO 

(A* Af,f) 


(75) 


(76i 


where the symbol + "+" denotes the transpose of the adjoint operator of the vector operator A- The 
product of the A A is a scalar operator and it is given by 


A A - A A. + A , A 2 +... + A „ A „ 


(77) 


Because the Af JL are the elements of the Hilbert space G, therefore the extreme value problem 
(75) is the same as the problem defined in (ibid, chapter 4). The solution of eq. (75) is then 


(A 'A + cf- Lm L m ) fa — A g 


(78) 


If the inverse of the operator 


A* A + u 2 L* U 

exists, then eq. (78j can be written in the form 
fa = (A + A + a 2 Lm U,)" 1 A + • g 

By using eq. (77) and the definition of the vector operator eq. (78) can be written as 


I 

i -1 


A , + 


2 . 
a L r 


L r 


/ 



i -1 


(79) 


(SO) 
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Applying the Fourier transformation to eq. (80) and using eqs. (22) and (64), we obtain the 
spectrum of the regularized function f a : 


ii * 

!♦. 


Cl). 


6 , 


i=l 


co f - 

1 a 


II * J * 

XM, + Ct $lA. 

i=l 


(81) 


In order to get the f a in space domain, the operator equation (80) has to be solved. In the 
frequency domain this problem becomes much easier: the regularized function f a can be obtained 
by taking the inverse Fourier transformation of eq. (81). 

The estimation error covariance is the same as eq. (71): 

C ce= C ff- Cfj-C ffa +Cf o f o (g2) 

^here the superscript "R" distinguishes the estimation error covariance between the heterogenous 
data and the homogeneous data cases. 

Notice that 


F !Cf,' = R 


f! 


F{C fJ }=O 0 R f . 


F;C. ta j = O 0 Rn 


f ' } C tV,( ~ °(.°(l R ff + °e • 


wnere 



n • 

l-l 


n • 2 • 

Z°, 0, + a OlJDl 

i=i 


(83) 


(84) 
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II • 

X°, °, R . 


1=1 




V 0 , o, + a 0 L J? U 

W --1 


then we nave 


R cc = U -O 0 ) Rff + 0 C 


( 85 ) 


( 86 ) 


obviously the Ree is a function of the parameter a and the power spectral density function Ol^- 


5. Smoothing 

In practice there are always some kinds of smoothing being used in the numerical 
computations. For instance, the use of the mean values of the data is a smoothing. A smoothing 
procedure can be used for solving the improperly posed problem. If we know the frequency 
composition of the function f and of the measurement error, we can design a smoothing operator 
i filter) to filter out the effect of errors. 

Here we introduce one smoothing operator which has the spectrum 

0 ,. = — --. ex > 0, X. > 0 

A 

1 + (l0) (87) 

where 0 S is the spectrum of the smoothing operator S, and a, X are parameters which can be 
chosen; 0 ) = (u- +■ v 2 ) 1 ^, u, v, are frequency variables. 0 S is a low pass filter because it filters out 
rhe high frequencies and lets the low frequencies pass through. 

For the improperly posed problem we have the solution 

f s = S A' 1 g (88) 

In the spectral domain eq. (88) can be written as 


0) f = 0 S 0 A o i 


(89) 


The smoothing procedure can also be applied to the heterogeneous data. If the data are 
errorless, the best estimation operator H, was given by eq. (53) and the spectrum of the solution is 
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( 90 ) 


n • 

lo, 

1=1 

n • 

y o,o. 


i 


in reality the errors are always included in the data. Except for the systematic error the 
measurement error are often modeled by random error and such error effects mostly the high and 
very high frequencies of the data Therefore a low-pass filter can be used to decrease such effect. 
Another reason of the use of the low-pass filter is that the high frequencies (nearby Nyquist 
frequency) must be minimized in the numerical computation. Because such frequencies are mostly 
distorted by the measurement error, sampling error, truncation error, etc. 

After die smoothing of eq (90) we obtain a solution: 


^>0 0 ) 
X- v i < 

! 1 


(!).- - 0. (!), 0 


n • 

\ 


0 0 


(91) 


If the ^mo*'trv.ng procedure ;s chosen properly, the effect of the errors can also be minimized 

a .rr.'i>th solution can tv; *>tained from an improperly posed problem. 

I"' ■: moot.’ur g error is given by 


(92) 


r ■ anatwe 


IS 


\1 ! e 'P| e' iQ)/ 


( ff - ( .. - C j f +■ C t i 


(93) 


j'lrtg e.|v (i 1 ), iSSi and (92 1 we get the power spectral density function of the 

■ 'trim; er. >r 


H,si -- 1 1 O j 9,; 4 - '0,0 A I R nn 


! >r - he reUToceneoiis data the estimation error e s is given by 


(94) 
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where R„ is the power spectral density function of the .neastirement error e;. 


( 95 ) 
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(i Relationship Between the I^ast Square Collocation, Regularization and Smoothing 

The relationship between the least squares collocation and the regularization method was 
>hown in (Rtimmel, et a!.. 1979) Both of the methods are identical under seme restrictions. 
Basically, three methods are the same and have the same property: they filter out the high 
frequencies in the solutions and nuke the solution stable and smooth. We will show that they are 
identical it they satisfy a few conditions. 

The easiest way to show the relationship between the three methods is to inspect them in the 
spectral domtin. \Ve consider only the homogeneous data case. Comparing eq. (32) with (69), 
we find that On,, and Cq become identical, if the following equation holds 


U 0 I-„°L„ = R nr/R|! 
Rewrite eq 1 69) in the form 


(96) 


0 M = 0 A \ 1 i- a 0, Oi 0. v 0 A t 
and compare (97) with (89) and using eq. (87), we get the equation: 


(97) 


u Oj 0, |0 A 0 A 


= a,cu 


(98) 


If the regularization parameter a and the spectrum of the operator L™ satisfy eq. (97), then the 
regularization method is ide uical with the smoothing method which has a spectrum similar io eq. 

i 8 7 , 

Rutting eqs. (96) and (98) together, we get the identical condition of three methods: 


■1 0 0: -- K R If - 0 A 0 A 0) 


(99; 


i: the regularization parameter u, Miiociimg parameters a s , X and the regularization operator 0i. m 
a:; };, eq t‘>‘* j. then the three methods are identical. 

1 he .:udy has -h..wn the relationship between the three methods. They pose the same property - 
.cir.g <>r filtering out the high frequencies in the solution. But they are different procedures 
I he;, are identical only when all of them fulfill the condition (99). 

for the regularization and smoothing method we can choose the parameters and the 
r' c u inzaiion operator 1 m or different smoothing operator to get snxxith solutions. Therefore they 
more flexible than least squares collocation for the solving of improperly posed problems. 










7. Numerical Test 


In this section we take numerical tests. The goal of this numerical simulation is: To have an 
idea about the use of the Taylor series to get the gravity disturbance on the earth’s surface by 
processing the aerial gravity-gradient data. The formulas above derived are used to determine the 
derivatives of the disturbing potential, such as Tz. T^, T zzz at the mean level. We want to know 
how good the methods are and have a view about the magnitude of the terms, such as Ah T zz , 
l/2(Ah) 2 Tzzz. 

The test area that was chosen has the geographic latitude 32°<(p<35° and the geographic 
longitude 257°<3.<260°. The 4 km x 4 km free-air anomaly for the United States (Rapp, et al.. 
1988) was used as the original data. 

In the computation the point mass model was used. We assumed that there is a point mass 
laver embedded at a depth of 8 km below sea level (i.e. the geoid). The relationship between the 
point layer and the disturbing potential T is given by: 

N m M 

T (x, y, d) = - XX -- - YJi 

,=lM [l*-^ 2 ^y-yj) 2+d2 l (ioo) 

where M,j is the product of the point mass at point (x,, yj) times the gravitational constant G, and d 
is the height of the computed point; N, M are the grid numbers of the area along x and v directions, 
the minus sign in eq. (1 (X)) is for convenience. The geometry of the point mass model is shown in 

Figure 2. 

z 



Figure 2. Geometry of the Point Mass Model (zq = 4 km, hm = 1.5 km) 
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The computational process is described generally as follows: 


1. Read the free-air anomaly from the 4 km x 4 km data file (Rapp, et al., 1988) and assume 
the data are given on the sea level. Furthermore, the free-air anomaly is assumed to be the 
gravity disturbance T z , and the formula (102) given in paragraph 7.3 is used to get the Mjj. d 
was chosen as 8 km; 

2. computing T„, i, j = 1, 2, 3 on the flight level by using formulas (103), (104), (105 
(107), (108), and (1101 given in paragraph 7.2; 

3. corrupting T,j by adding random errors with error variances c 2 = 1, 4, 25 Eorvos 2 and 
mean value equal zero; 

4 processing the simulated gradient data to get T„ T^, T 1ZZ on the mean elevation level; 

5. comparing the computed i",, f lz , ^\ 1Z with the "true" value which were computed directly 
from the point mass model. 

In the following we give the description of the data, the formulas used and some considerations 
about the numerical computations. 


7.1 Data Used 

The gravity anomaly in 4 km x 4 km grid point values for the United Slates was used. In 
order ;o get higher frequencies in the solution, the data was interpolated in 2 km x 2 km gnd 
interval by using the bicubic spline function. All computations were based on the data in 2 km x 2 
Mr gnd values. The use of 2 km x 2 km gnd interval is based on the following considerations. 

1 Using a small grid interval can decrease the aliasing effect in the Fourier transformation, even 
though we do not get more information by interpolating the 4 km x 4 km gravity anomaly into 2 
km x 2 km gnd point values; 

2. The evaluation of the formulas, such as given in the next paragraph, can be more accurate by 
using 2 km x 2 km gnd than using 4 km x 4 km grid. Using the smaller gi id can increase the 
computation accuracy (cf. Tziavos. ct al., 1988). 

The statistics of gravity anomaly Ag in this area is given in Table 1. 

Table 1. Statistics of the Gravity Anomaly Ag in the Test Area. 

(mgal) 



j RMS value 



-8.X2 

i 24.48 

78.21 

-76.17 


The contour map of the gravity anomaly in the test area is shown in Figure 3. 
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7.2 Formulas Used 


Before the numerical tests were taken, the formulas which would be used are written in the 
follow ing. The spectra of the differential operators used occur in the airborne gravity gradiometry 
are given in Table 2 (cf. Vassiliou, 1986). The definition of the spectra of the operators are given 

in eq. (22). 

Table 2. Spectra of the Differential Operators 


Operators 

Spectra 

a 


dx 

j 2 n u 

a 


dy 

j 2 k v 

a 


az 

- 2 Jt 6) 

n 

a 

k 1 p 

dx dy dz 

(j 2 n u) (j 2 re v) (- 2 it to) 


Here we have k + 1 + p = n; k, l, p = 0, 1,2,... 

It is easy to find the spectrum of the upward continuation operator U defined by eq. (1): 

0 =e _:W “ 

°u e non 

The relationship between the disturbing potential T and the mass point M,, is given by eq. 
1 100). The derivatives of the disturbing potential T,, T,,, i, j = I, 2, 3 can be aerived from eq. 
(100) (cf. Vassiliou, 1986): 


M N 




i-i j=i r 


T„=IX 

,=i j=i 


M N , ,2 , .2 .,2 


j*-*,) +ly-yj - 2d' 

,5/2 


M, 


V' v' 3 lx x J d 

t, = - y y —s- m,, 

/x ,5/2 1J 


'-'I C l 


( 102 ) 


(103) 


(104) 
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V- £ 3 (y-yjd 

i = 1 j= 1 1 * 


M N ¥ v 

t =y y —- m, 

X. i—d d—d -111 1 . 

1=1 J=1 ) 


M N 


i = l i=l 


IY'-Vj) + d" - 2 (x-xi) 


,5/2 


■M, 


T 


*y 


-II 

i=i j=i 


3 (x-x,) (y-yj 


,5/2 


M, 


( 105 ) 


(106) 


1 r\-i \ 
V *) 


(108) 


M N 

■-II 

i = l )=1 


y-yj 

,3/2 


M; 


(109) 


vy i—t i—d $12 ‘1 

>= 1 j=i i (HO) 

where I = [(x-x,) 2 + (y-yj) 2 + d 2 | l ' f2 . 

Taking the Fourier transformation of the formulas (102) - (110), and using Table 2, we have: 


F |T) = F {M.jJ 



F ( T Z J = F f M, j } (2 7t e 2n0)d ) 

F{T / ,) = F(M 1J ){-(2 I t) 2 coe 2 ’ tox1 / 
F (T,x) = F {j ( 2 k )* ue 2n<od } 

F jT/y) = F j M,j) (j (2rt) 2 ve 2na>d / 
p l T i = F i M 

' X( 1 1J M J oi / 


( 111 ) 

( 112 ) 

(113) 

(114) 

( 115 ) 

(116) 
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(117) 

(118) 

(119) 

( 120 ) 


Here we should not confuse the subscript j with the imaginary number j = V-T. 

In the numerical computations eqs. (Ill) - (120) were discretizated and were evaluated by 
using the fast Fourier transformation (FFT). The first and second derivatives of the disturbing 
potential were computed at the flight level and mean elevation level by using eqs. (111) - (120). 


After the computations of the second derivatives of the disturbing potential T,j at the flight 
level, the normal distribution of the random noise with the variance o 2 = 1,4, 25 Eotvos 2 was 
added to the computed T,j. In processing the corrupted data to determine the derivatives of the 
disturbing potential Tj, T^, l uz on the mean elevation level, the regularization method was used. 


It was assumed that the disturbing potential T and its partial derivatives to z up to third order 
are smooth at the mean elevation level. For the word ''smooth'' we understand it under such 
meaning: The disturbing potential T has continuous partial derivatives up to third order. In fact, 
the disturbing potential T has continuous partial derivatives up :o all orders above the earth's 
surface. Because the analytical downward continuation is used, we constrained the disturbing 
potential has continuous partial derivatives up to 3rd order above the mean elevation level. We 
took the regularization operator in eq. (75) as 





( 121 ) 


where U is the upward continuation operator. 
The spectrum of the operator L m is given by 

0 L - - (2710)) Oy, 


( 122 ) 
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Six components of the second derivatives were used to determine Tj, i = 1, 2, 3 on the mean 
elevation level. The observation equations are 


L = A 

* 

T. + £ 

* X 

(123) 

1 = A 

T v + £ 


Y 

Y y 

(124) 

L = A r 


(125) 


where j_ is the observation vector which has fhe components Ij, i = 1,2, ... 6, the measured second 
derivatives of the disturbing potential, Aj, Ay and are vector operators; £ is the vector of the 
measurement error. As a specific example eq. (123) is written in the form: 


yz 

l 




yy 


A 2 

A 3 


Tx + 


(123’) 


where tJ z , Ty Z , ... are the measured gradients of the gravity disturbance; ei, E2.••• are the 
corresponding measurement error. 

The best estimates of the T x , T y and T z in eq. (123), (124) and (125) are given by: 

T x (x, y- h m) = F ' (k|u, v)j 2nu} 


T y (x, y, h m | = F ' (k|u. v)j 2 kv} 
T/( x . y. h rT J = F 1 {k (u, v)(- 2 ko))| 


(126) 

(127) 

(128) 


where the symbol " A " denotes the estimate of the function, and k (u, v) is given by: 
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• 1 6 • 

$u W g, 

1=1 

k(u, v) =- 

6 » 2 I 

4>i+ a (2itco) 
i=i 


( 129 ) 


Here we have used eqs. (81) and (122). The co g j is the Fourier transform of the observations T^. 
From eq. (126) - (128) one can see that the k (u, v) is nothing but the spectrum of the disturbing 
potential T determined by using the gradient data Ty. If the Tjj is ordered as in (123'), then we 
have: 

2 2 
9l = — j (27c) uco, <{) 2 = — j (27c) VO), 


0 3 = (2ti) CO , 0 4 = -(2ji) U 2 , 



Insert (130) into (129), we get 


k (u, v) = 



(130) 


(131) 


In the numerical computations the regularization parameter a has taken the value 0.07. 

In the same manner we have the formulas for the second and third derivatives on the mean 

elevation level: 




T„.(*.y.hm) = F ' 

'f k , 

|u, v)(2 tt) 

(-juco) / 

(132) 

T,,!.x.y.h n j = F ' 

!ki 

(u, v}(2it) 

f jvco) / 

(133) 

~ .] 

/ 

C 

< 

g 

to 

A 


TJx.y.h^F 

1M 

CO / 

(134) 


,/ 

|u, v) (2 Tt) 

W 2\\ 


T„ / (x.y.h n J = F- 


'JUCO ff 

(135) 
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( 136 ) 


T vw (x. y, h ni ) = F' 1 {k(u, v)(2jc) (jvco )} 
T 77i (x, y. h ri J= F 1 \k(u, v)(2k) ( co )} 


(137) 


The formulas (126) - (137) were discretizated and evaluated by using the fast Fourier 
transformation. 


If we use the smoothing method, according to the equation (98), the parameter a s , \ are 
chosen 


[ 2 

a s = - a = 100, X = 2 


The spectrum of the smoothing operator is 


(138) 


1 



1 + acu 1 + 100 a) 


The least squares collocation method can also be used. But the power spectral density 
function of the disturbing potential and of the measurement error have to be properly chosen. 

The regularization error was not computed by using the formulas like (86), because the "true" 
value could be computed directly. The difference between the "true" values and the computed 
values were computed and the results are shown in the following section. 

If we have no true" values, as it should be in practice we can choose the spectral density 
function of the disturbing potential T and of the error E, the estimation error (regularization error, 
smoothing error) can be estimated by using the above derived formulas. 


7.3 Considerations of the Singularity and Recovery of the Spectra of the Gravity Disturbance 
From Gradient Data 

In the above derived formulas we have the singularity problem. The formulas, e.g. eq. (131). 
are not defined at the origin. Such a problem can be catalogued by the singular integral problem. 
A theoretical study of singular integrals and integral equations is given in Miklin (1965). The 
study of the singular integrals in the physical geodesy can be found in (Siinkel, 1977; Wang. 
1986). 

Here we consider the singularity problem in the spectra domain. As a common example, we 
consider the Stokes integral. In the planar approximation we have: 


N = 


1 


2 try 



Ag 


dx dv 


(140) 
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where N is the geoid undulation, y. and Ag are the normal gravity and gravity anomaly, 
respectively. Applying the Fourier transformation to eq. (140) we get 


o) N= -i_to g 

(141) 

where con, Wg are the Fourier transformation of the geoid undulation and the gravity anomaly, 
respectively. 

We have two methods to compute eq. (140). One is taking the discrete Fourier transformation 
to eq. (140). The kernel T 1 has a singularity at the point (x = 0, y = 0). The treatment of the 
singularity can be found in (Heiskanen and Moritz, p. 121; Schwarz, et al., 1989). 

Eq. (141) can also be used to compute the geoid undulation. The only question is, what kind 
of value eq. (141) should be taken at the point (u = 0. v = 0)? It will be shown that the choice of 
the value coy; (0, 0) effects on the geoid a constant bias. 

Let 

co N (0, 0) = (3 (0^ (O, O) * 0 


where P is an arbitrary constant. Notice that 


(o (0.0) = 


// 

II 


oo 


2nj (0 **(} y) 
e 


Ag (x, v) dx dv 


Ag (x , y) dx dy 


(143) 


If the integral in eq. (143) is limited in a local area, as the case in practice, cog (0, 0) is not always 
equal to zero. Assume that eq. (141) is discretized and inversed by using the discrete Fourier 
transformation, then the o)\ (0, 0) has the contribution to the geoid. 
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n - 0 n 0 


w 8 (0.0)< 


t m 0 
Zk > TT' 


n ()\ 
N 


= B oj (O, 0 ) 


(144) 


where 5N is the change of the geoid due to the different choices of the con; (0, 0). 5N is a constant 
every where. If we set cu g (0, 0) = 0, then 5N = 0. Therefore the geoid undulation from eq. (141) 
may have a constant difference with the geoid undulation direct from the Stoke's integral (140). In 
order to r move this bias, other data, such as reference fields, should he used. 
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Another question is: Can we recover all spectra of the gravity disturbance from the gradient 
data 9 The answer is no. For example, the relationship between T z and T zx is 


0) T = j2uco T 

1 a 4 i 


(145) 


If we want to recover from we have 


. 1 

CO'j' — J 03-r* 

1 z 1 i * 

2lZu (146) 

At the line u = 0, the spectra (oqy are not defined. A common way in the numerical computation is 
to set 

0)j 7 (0, v) - 0, v e (- «, + oo). (147) 

Butin reality eq. (147) is not correct, since the spectra cur z (0, v) should not be equal to zero. The 
conclusion is. The spectrum of T z cannot be recovered from the spectra T zx entirely. 

In the numencal computations, we have to take the assumption (147). Therefore it is not good 
enough using T xz (including T yz ) to determine T, because we cannot recover all frequencies of T z 
from T /x or T zy . 

A better recovery of the spectrum of T z can be achieved from T^- 


OJ 


T,~ 


- co , 


2ttO) 


(148) 


This function is not defined at the origin (0, 0). If we assume the mean value T z is equal to zero, 
based on eq. (143) we have: 


oi T (0,0) = 0 

1 t 


(149) 


and we can get all frequencies of T z from T u . If the copy is not equal to zero, we still have to take 
(149) in the numencal computation. But the information on the bias of T z has to be provided by 
other data, such as the point disturbance component values on the ground (Jekeli, 1986). 

The best way to recover the frequencies of the gravity disturbance is the use of all components 
of the second derivatives of the disturbing potential. From eqs. (128) and (131) we can see, 
frequency of T, from all data is not defined only at the point (0, 0). We define the frequency of T 7 
at this point as equal to zero, and the frequency copy (0, 0) is obtained from tie point disturbance 
component values or another data. 

Summary: We can recover all frequencies of the gravity disturbance by processing the aerial 
gravity-gradient data except the mean value. The mean value of the gravity disturbance must be 
provided by tie point values or another data. 
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"4 Rc>uits 

!n ihe pax ceding action % ■■•■■' dew ribed the data, the formulas used. In this section we give 
: .c computational row is V.' • di : . aiivc' of the disturbing potential were determined on the mean 
elc o:.’n level a.i.i coo,pared wnr. the true" values which were computed from the point mass 


I able 3 we give the statistics of the differences (errors) between the computed and the 
true values for the partial derivatives of the disturbing potential to z. In the numerical 
comp.nation eqs i ;28i, < 134) and < I -7 ; were used. The regularization function was chosen as 

described in section 7 2 


Fable 3. Statistics of the Ddfercnccs Between the Computed and the "True" Values 


j 

Noise 'eve! <h) 

mean 

RMS 

max 

min 


o' - ! 

- 0.29 

0.52 

2.57 

- 3.72 

1 ( ^, . 

a ~ 2 

0.29 

0.60 

2.78 

- 3.94 

} 

a---5 

- 0.29 

0.99 

3.67 

- 4.76 


c - 1 

- 0.03 

0.26 

1.47 

- 1.64 

i : -c 

o = 2 

- 0.03 

0.47 

2.13 

- 1.95 

1 ^ « ri y . 

a = 5 

- 0.03 

1.14 

4.35 

- 4.73 

j 

a= ! 

0 84x10^ 

0.40 

1.71 

- 1.69 

i I/v 

o = 2 

- 0.87x10-3 

0.77 

3.05 

- 1.95 

! ’V.gaiT m-, 

n = 5 

- 0.96xl0- 3 

-L21L 

7.78 

- 6.97 


From Table 3 we can see that the random noises have no significant effects on T z . In the 
..ornpc:,irons the regularization was used and the effects of the noise have been mostly removed. 
T: .s not the s«rne as expected iJekeli, 1^87) The reasons may be: 

1 We u-.-d the gravity anomaly in a not very rough (gravity anomaly wise) area. The spectra 
: •n'.v.ed mostly by lower frequencies. The power spectral density function of the data is 

•V": n Figure 4 The amplitude of the high frequencies of the data is small and should have no 
'.:i: tic art c; -rtr:button to i 'X u and T m . The high frequencies of the measurement noise and the 
:..c:i ’reqrrrcie- >./ the data were filtered out and the results were not changed significantly from 
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Figure 4. Power Spectral Density of the Gravity Anomaly Ag ; n Test Area 

2. An important thing in the numerical simulation is the choice of the grid inteval. The 
gradient data used were regularly distributed in a 2 km x 2 km grid interval. The small gnd interval 
minimized the aliasing effect. A study about the computation accuracy with the grid spacing, flight 
altitude can be found in (Tziavos, et ah, 1988). Even though where they were talking about the 
computation of the topographic effect on the gradient data, we get some idea about the relationship 
between the computation accuracy and the ratio of the grid spacing/flight altitude. In our 
computation the ratio of grid spacing/flight altitude was 2/2.5 and it met the computation accuracy. 


3. The rauo of signal/noise is defined by the RMS value of signal/RMS value of noise. In the 
numerical computations this ratio was reasonable. For example, for the gradient T u we have the 
signal/noise ratio equal 8, when the random noise has variance a = 1 Eotvos. A very interesting 
phenomenon is, even if a = 5 Eotvos when the ratio of the signal/noise equal 1.6, the results for T z 
are still good. The reason may be: the random errors pose high frequencies and they are filtered 
out by using 'he regularization method; the second reason is that all components of the gravity 
gradient have been used and it improved the results. 

Of course we cannot expect to obtain good results of T zzz by processing ihe aerial gradient 
data. The high derivatives of the disturbing potential T zzz is more sensitive to the measurement 

-rr< ir 

In the following we determine the derivatives of the disturbing potential up to third order on 
mean elevation level, and the Taylor's series are used :o get the gravity disturbance on the 
earth's surface The difference between the "true” values and the computed values are given in the 
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following tables. The second correction term of eq. (9) is also included and the results are 
presented in the tables. The results w tthout the second correction term of (9) are also shown in the 
following tables. 


Table 4. Statistics of the Differences of 


t 2 - 


T z+ Ah T u -4(AhrT ZZ7 J 


unit mgal 


noise level (E) 

mean 

RMS 

max 

IEZQH 

0 = 1 

- 0.06 


7.25 

- 7.55 

0 = 2 

- 0.06 


8.52 

- 9.08 

0 = 5 

- 0.06 

WBM 

18.83 

- 17.54 


Table 5. Statistics of the Difference of 

vfh+Ah 


unit mgal 


noise level (E -1 

mean 

RMS 

max 

■sm 

a = 1 

- 0.06 

0.91 

7.89 


0 = 2 

- 0.06 

1.26 

8.51 

- 9.43 

0 = 5 

- 0.06 

2.61 

10.39 

- 12.04 


Comparing Table 4 with 5 we find: If the accuracy of the gravity gradient data is poor, e.g., 0 
= 5 E, the term 0.5 (Ah) 2 T JXZ is corrupted by the errors and it is not usable (compare the results in 
Table 4 and 5 when 0 = 5 E.). But this term still gives some contributions to the big values of T z 
: compare the results in Table 4 and 5 when o = 1 E.) while the RMS values becomes bigger. It is 
expected that this term gives significant contributions to the big absolute values of T z if the gradient 
data are accurate and in good distribution. 

The map of the contour line for the estimated *t 7 was drawn in Figure 5. Comparing Figure 3 
with 5. one can find the T z is smoother and smaller (absolute magnitude). Figure 6 shows the 
contribution of the correction term Ah 't' zz , in test area. A few significant corrections are in the 

rough gravity anomaly area. The contribution of the correction term l/2(Ah) 2 i* zzz is shown in 
Figure 7. The correction is small and rough. This correction can become very small after some 
kinds of smoothing, e.g., the results are averaged into mean block values. Figure 8 gives the 
computed T z in the test area. In comparison with Figure 3 one can say that the recovery of the 
gravity disturbance by processing the gradient data are successful. Figure 9 gives the difference 
between the ' true'' values and the computed values. Although the difference (error) can reach 5-7 
mgal at some points, but it can become smaller when the results are averaged into mean blc k 
values. 
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Figure 6. Contour Map of Ah *t /y in Test Area 
Contour Interval = 1.5 mgal 
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Figure 9. Contour Map of the Difference T, - ( + Ah ^+ 1/2 (Ah) 2 ^tzu) 

Contour Interval = 2 mgal 
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Similarly, the computation at results for T x , T y are given in Table 6-11. 

Table 6. Statistics of Differences Between the Results and the "True Value" 
unit mgal 



Table 7. Statistics of Differences of 



unit mgal 



Table 8. Statistics of the Differences of 
T x -(f x + Ah T zx ) 


unit mgal 


noise level (E) 

mean 

RMS 

max 

mssm 

o = 1 

0.G1 

0.63 

6.48 

- 5.67 

n 

II 

O 

0.01 

0.88 

6.56 

- 5.77 

0 = 5 

0.03 

1.84 

7.41 

- 7.44 
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Table 9. Statistics of Differences Between the Results and the "True Value" 
unit mgal 



Noise level (E) 

mean 

RMS 

max 

min 

T y 

O = 1 


0.31 

1.97 

- 2.89 

(mgal) 

0 = 2 


0.38 

2.19 

- 2.94 

0 = 5 


0.67 

2.85 

- 3.35 

■mu 

o= 1 

- 0.90x10 

0.18 

ifHEl 

- 0.96 

1 1 

0 = 2 

- 0.87x10-3 

0.33 

■ 

- 1.30 


0 = 5 

- 0.79x10-3 

0.81 

EBJ 

- 3.18 


0 = 1 

-0.17x10-”? 

‘0.28 

1 06 

- 1.15 

Tyzz 

(meal/km-) 

o = 2 

0=5 

- 0.21x10-3 

- 0.32x10-3 

0.54 

1.35 

2.09 

5 18 

- 2.14 

- 5.37 


Table 10. Statistics of Differences of 


T - 

y 


/ i 2 \ 

T y + AhT„ + 'l-\h) T,„j 


* 


unit mgal 


noise level (E) 

mean 

RMS 

max 

min 

0 = 1 

0.11 

0.79 

4.52 

- 5.86 

0 = 2 

0.11 

1.38 

5.55 

- 6.07 

0 = 3 

0.12 

3.27 

12.54 

-12.81 


Table 11. Statistics of Differences of 
V(f v + AhT y/ ) 


unit mgal 


noise level (Ej 

mean 

RMS 

max 

min 

O = 1 

0.11 

0.65 

5.02 

- 7.21 

0 = 2 

0.12 

0.89 

5.42 

- 7.25 

0 = 5 

0.13 

1.85 

6.99 

- 8.19 


From Tables 4-11 wc come to the following conclusion: The gravity disturbance can be 
determined on the earth's surface with demanded accuracy. With 1 Eotvos error in the gradient 
data, the components of the gravity disturbance T x , T y and T z are determined with an accuracy in 1 
mgal. Even though the measurement error is 5 E., the gravity disturbance can be determined with 
an accuracy of 3 mgal on the topographic surface. 

Here we need to point out that the measurement error model is assumed a normal distributed 
random noise. Although this assumption is not entirely realistic it gives the main property of the 
effects of the error in the processing of the gradient data. It is expected that if the spectra! 
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distributions of the gradient data and the measurement error are known, one can use the methods 
proposed in this study to minimize the effect of the measurement error and to get a stable and 
reasonable solution. 

Normally, for a convergent senes, the more correction terms that are taken, the results are 
determined more accurately. But we do not know whether the Taylor's series is convergent on the 
mean elevation level. From Tables 4-11 we can see that the correction term l/2(Ah)*T zzz may 
have some contribution to the large values of the gravity disturbance when the measurement error 
is small. But it supplies wrong information when the data accuracy is poor. Therefore we suggest 
if the area is very' rough and the data are accurate enough, e.g. measurement error is lower than 1 
E., then the correction term l/2(Ah) 2 T zzz , etc., could be taken into account; in a flat area this term 
is very small and can be neglected; if the accuracy of the gradient data are poor, the correction term 
l/2(Ah)- T zjz, etc., could be wrong and it is not of benefit to the results. 


8. Conclusion 

After many years of development the airborne gravity gradiometer survey system is coming in 
to practical application. Recently, a test flight was taken in the Texas-Oklahoma area which is 
characterized by a very smooth topography. Although this was the first time the gravity 
gradiometer survey system was flown, and only a fraction of the total runs yielded good gradient 
data, the components of the gravity disturbance were determined on the ground with the accuracy 
of 2 to 3 mgal in 5' x 5' mean anomalies. 

In the future the test will be carried out in the rough mountain area and the topographic effect 
has to be taken into account. This problem has been considered by many authors. Tziavos et. al, 
have developed an estimation algorithms for the computation of the effect of the mass above the 
ellipsoid (Tziavos, et al., 1988). If we subtract the contribution of the topography from the 
gradient data, it means the mass of the earth is adjusted by removing the mass above the ellipsoid. 
This is not correct in some cases, e.g. the determination of the geoid. This report did not discuss 
this problem. Our goal was the determination of the gravity disturbance on the earth's surface. 
We assumed that the disturbing potential and its derivatives can be analytically downward 
continued to a mean level - in the report the mean elevation level was chosen, then the Taylor series 
was used. The gravity disturbance was given by (cf. Figure 1): 

T Z (QJ = T z (p') + Ah T JP') + i(Ah)' T„ 7 (p'} 

Obviously this analytical downward continuation problem is an improperly posed problem. The 
solution (Q) may pose serious numerical difficult and not be stable. 

For an improperly posed problem there are three methods that can be used: the least squares 
collocation, regularization and smoothing. How can they be used in the gravity gradiometry was 
studied in Section 3, 4 and 5. The studies indicate that the three methods are essentially the same. 
They filter out the high frequencies of the data and make the results smooth and stable. In 
comparison to the least squares collocation the methods of the regularization and smoothing are 
more flexible. 

The regularization was used for a simulated computation. The numerical computation shows 
that this method is qualified for the analytical continuation of the derivatives of the disturbing 
potential, such as T z , T zz , T zzz to the mean elevation level. When the accuracy of the gradient data 
was 1 Edtvos, the gravity disturbance was recovered with the accuracy of 1 mgal. If the accuracy 
of the gradient data is poor, but the ratio of the signal/noise is still high enough, the recovery of the 
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gravity disturbance can achieve a reasonable accuracy. One example described in the report is 
when the gravity disturbance was determined with an accuracy of 3 mgal when the accuracy of the 
gravity gradient data was 5 Eotvos. 

The derivatives of the disturbing potential T,. ; , can be obtained by processing the aerial 
gradient data but it is corrupted by the measurement errors. It is still a difficult work to get the 
higher derivatives of the disturbing potential and sometimes it looks like it is impossible. 
Fortunately the high derivatives of the disturbing potential have most effects on the high 
frequencies of the disturbing potential which do not have significant contribution to the disturbing 
potential. 


For a flat area the first two term' eq o mls/ied our needs, but in n rough mountain area 
the last term in (9) has to be considered if the accuracy of the gravity gradient data is high. If the 
accuracy of the gravity gradient data is poor, Ge iasi ut.m is not beneficial to the results. 

The numerical simulation war limited by different assumptions, such as the measurement error 
model, no position errors in the data, etc., but it ;cpresents the main property of the processing of 
the aerial gradient data and appro , ; . d world to a t .:eat extent. All the computations 

were completed by using the fast Four or transf onmtio'i. 
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